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Abstract Mixed finite element methods solve a PDE using two or more variables. The theory of Discrete Exterior 
Calculus explains why the degrees of freedom associated to the diflferent variables should be stored on both primal and 
dual domain meshes with a discrete Hodge star used to transfer information between the meshes. We show through 
analysis and examples that the choice of discrete Hodge star is essential to the numerical stability of the method. 
Additionally, we define interpolation functions and discrete Hodge stars on dual meshes which can be used to create 
previously unconsidered mixed methods. Examples from magnetostatics and Darcy flow are examined in detail. 



^ 1 Introduction 

o 

^SJ ■ The theory of Discrete Exterior Calculus (DEC) has provided a novel viewpoint for analyzing linear systems derived 

from finite element theory. We highlight three important conclusions of this theory; 

^ 1. Variables in a PDE should be discretized as degree of freedom arrays ("cochains") over a primal simplicial mesh 

or its dual mesh. 

2. A discrete Hodge star is used to transfer information between primal and dual meshes. 

3. Whitney elements provide stable finite elements for the primal mesh. 

Most numerical methods for PDEs over unstructured tetrahedral meshes discretize variables as cochains over the 
, primal mesh and build up linear systems from there. In this paper, we look at the alternative approach of discretizing 

variables over the dual mesh and design dual formulations of the linear systems based on DEC theory. This approach 
is especially valuable in the context of mixed finite element systems as they employ all the key ingredients of DEC 
• theory: both primal and dual cochains, a discrete Hodge star, and, typically, Whitney elements. 

^ Before turning to mixed systems, however, we look at a simpler example from electromagnetics illustrating the 

relevance and benefit of our technique. The example is inspired by He and Teixeira [T7]. Using a Discrete Exterior 
Calculus analysis of Maxwell's equations, one can derive a second order vector wave equation 

'sj- ■ M2lD)lE = OJ^MiE, (1) 

> , 

where E is the electric field intensity, discretized as a cochain on the primal mesh, tj is a coefficient, Di is a rectangular 
incidence matrix having entries of and ±1 only, and Mj, is a discrete Hodge star operator. 
■ The dual formulation of this physical phenomenon is an equation for the magnetic field intensity H, discretized 

as a cochain on the dual mesh: 

DiM^^dIh = Lo'^M^^'a. (2) 



(N 



Both systems ([T]) and ([2]) are computationally tractable if Mj. is a diagonal matrix which, by DEC theory, can be 
achieved when the primal and dual meshes are orthogonal. If orthogonality is not guaranteed, as is the case with 
barycentric dual meshes, Mj. is defined using Whitney elements and results in a sparse matrix. As a consequence, 
, system ((2|) then involves possibly full rank matrices and is thus significantly more computationally expensive to solve. 

5^ ' He and Teixeira [T7] reduce the rank of the M^T^ matrices by using a topological thresholding technique which requires 

an input parameter. 

Our approach skirts the problem of full rank inverses by introducing a novel definition of the M^^^ matrices free 
of parameters and guaranteed to produce a sparse matrix. The outline of the paper and summary of its contributions 
are as follows: 

— In Section O we briefly discuss prior work and fix relevant notation. 

— In Section[31 we use the Sibson coordinate functions to construct dual Whitney-like functions which define a novel 
sparse inverse discrete Hodge star (M^"°')~^. We show how the choice of discrete Hodge star requires certain 
geometric quality conditions of the primal and dual mesh elements. A specific example is given showing how our 
dual formulation of the problem can result in a better conditioned linear system than the primal formulations. 
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Fig. 1 Primal simplices are shown in black in the top row: cr" is a vertex, is an edge, and cr^ is a face. Their corresponding 
dual cells for n = 2 are shown in red on bottom: *a'^ is the barycenter of cr^, is an edge between barycenters, and Vrcr" is 
a planar polygon with barycenters as vertices. In three dimensions (n = 3), primal vertices have dual polytopes, primal edges 
have dual polygonal facets, primal faces have dual edges, and primal volumes have dual vertices. 



— In Section |4l we examine how our methodology applies to generic PDE problems as well as to some specific 
applications employing mixed finite element methods. We cast each into our common notational framework and 
show how to formulate equivalent dual formulations of the problem from a DEC-based analysis. The specific 
advantages of these dual formulations are analyzed, including an ability to compare and contrast calculations on 
a primal mesh with the analogous calculations on the dual mesh. 



2 Prior Work and Notation 

Our work is inspired primarily by the emergent theory of Discrete Exterior Calculus (DEC). DEC is an attempt to 
create from scratch a discrete theory of differential geometry and topology whose definitions and theorems mimic 
their continuous counterparts [19ll^ . A central conclusion of the theory is that degrees of freedom for finite elements 
should be assigned to mesh vertices, edges, faces or interiors according to the dimensionality of the variable being 
modeled. If these degrees of freedom have a natural geometric duality, as occurs for example between electric and 
magnetic fields, two meshes of the domain are necessary - a primal and dual mesh [TS] . This has given rise to DEC- 
based methods for solving problems of Darcy fiow [50], electromagnetism [T7] and elasticity [23, among others. As 
we will show, the 'bottom-up' approach of DEC clearly suggests alternative discretization methods less evident from 
such 'top-down' theories as finite element exterior calculus [2]. 

The main notational aspects of DEC are encapsulated by Figures [1] and [2] Figure [T] shows our notation for domain 
elements, i.e. primal fc-simplices cr*^ and their geometric dual n — fc-cells *cr"~'^ where n is the dimension of the 
domain. The dual domain mesh is defined by taking the circumcenters or barycenters of n-simplices and connecting 
them based on simplex adjacency in the usual manner. The measure of cr'^ (respectively ★cr""'^) is denoted |cr''| 
(respectively | ★(t"~''|), meaning length for k = 1, area for k = 2, and volume for k = S, with the convention that 
\a°\ = = 1. 

Figure [2] shows the various continuous and discrete spaces relevant to DEC theory for n = 3 and the operators 

between them. The vector space of fc-cochains, i.e. linear mappings from fc-simplices to R, is denoted C*^. The vector 



space of dual fc-cochains, i.e. linear mappings from fc-cells of the dual mesh to M, is denoted C . The Dj. matrix is 
the transpose of the {k + l)st boundary operator, i.e. it encodes element adjacency and orientation information with 
entries ±1. 
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Fig. 2 The combined DEC and deRham diagram for a contraetible domain in M^. The top row shows the deRham diagram 
with continuous Hodge star maps between function spaces. The middle and bottom rows show primal and dual cochain spaces, 
respectively, along with the discrete exterior derivative and discrete Hodge star maps. The I and V maps are interpolation 
(Whitney) and projection (deRham) maps. 



The interpolation map T}; converts a fc-cochain into a piecewise-defined fc-form whose global continuity in a 
distributional sense is indicated by Figure [2] (e.g. IiW G _ff(curl)). Define by 

2:fe(w):= w(a'^')W,.. (3) 

where is the Whitney function associated to simplex . These functions are described in Appendix A. The 

Whitney functions were first described in [23 and later recognized by Bossavit [6] and others as the correct gener- 
alization of edge and face elements needed for DEC theory. An extensive treatment of all of these spaces, functions, 
and operators is given in |14) . 

We now discuss the Hodge star * and its discretization as a square matrix M or M^^. As shown in Figure [2l the 
continuous Hodge star * maps between forms of complementary and orthogonal dimensions, i.e. * : yl*^ — > yl"^'^. For 
domains in as considered here, * is defined by the equations 

*1 = dxdydz, *dx = dydz, *dy = —dxdz, *dz = dxdy, ** — 1. 

For a more general definition of *, see [T]. 

A discrete Hodge star M maps not only between cochains of complementary dimensions but also between primal 
and dual meshes [18]. In this paper, we focus on the two definitions of a discrete Hodge star most relevant to DEC 
theory. The first is the diagonal discrete Hodge star defined by 

{M^^^%, - ^-^S^r (4) 

The definition of M^*°® fits nicely into DEC theory when the dual mesh is defined by taking circumcenters of the 
primal simplices, thus producing orthogonal meshes [9]. In practice, however, it is often desirable to use barycenters 
to define the dual mesh as this guarantees that will intersect -ka^ in the ambient space. A correction factor for 
this change is given by Auchmann and Kurz [3]. 

The more widely used approach for barycentric dual meshes employs Whitney interpolants in the definition of 
the discrete Hodge star: 

(Mf^^*),, (W,. , W^k ) = 1^ W,. ■ W^k (5) 

The inner product here is the standard integration of scalar or vector valued functions over the domain. Dodziuk [11] 
originally proposed the definition of M^'"* but it has been called the Galerkin Hodge [7] for its relation to finite 
element methods. Bell [S] has implemented linear solvers in a DEC context using MJ^''** for various k. 

Many other discrete Hodge stars appear in the literature, including the combinatorial discrete Hodge star of 
Wardetzsky and Wilson [25ll28j and the metrized chain Hodge star of DiCarlo et al. [10]. To our knowledge, no 
authors have defined a discrete Hodge star using dual interpolatory functions as we propose in this work. 
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3 Dual Whitney Interpolants and Dual Discrete Hodge Stars 

It is evident from tlie DEC-deRham diagram in Figure [5] that tlie direct interpolation of degrees of freedom on a dual 

mesh is not available in the common theory. Further, we have seen from the discussion in Section[5]that the definition 

of (Mfc)~^ has only been implied from definitions of Mj,. In this section, we define a set of interpolation functions X 

analogous to the Whitney functions and use them to provide an explicit definition of a dual discrete Hodge star. 



Define the dual Whitney interpolant of a dual fc-cochain w G C to be 

Tfe(w):= W(*a"-'^)W,,„-. (6) 

where W^g-n-k is a dual Whitney function associated to the fc-cell ★(7"-'= in the dual mesh. These functions are defined 
using a generalization of barycentric coordinates known as Sibson functions [22] , also called the natural neighbor or 
natural element coordinates [23 . Figure [3] summarizes the definition. 




Fig. 3 Geometric calculation of a Sibson coordinate, d is the area of the Voronoi region associated to vertex inside T. 
Z){x) is the area of the Voronoi region associated to x if it is added to the vertex list. The quantity -D{x) n d is exactly -D(x) 
if X = Vi and decays to zero as x moves away from , with value identically zero at all vertices besides v,; . The bottom right 
figure shows how the level sets of the Sibson coordinate associated to sit inside a single polygon. More figures can be found 
in Milbradt and Pick [21] 



Definition 1 Let x be a point inside a polyhedral cell T of the dual mesh. Let P denote the set of vertices {vj} and 
define 

P' = PU{x} = {vi,...,VAr,x}. 
Denote the Voronoi cell associated to a point p in a pointset Q by 

Vq(p) ~{yeT : |y - p| < |y - q| , Vq G Q \ {p}} . 

Note that these Voronoi cells have been restricted to T and are thus always of finite size. Fix the notation 

d := |Vp(vi)| = |{y G r : ly - Vil < |y - Vj | , Vj / i}\ 

= area of cell for in Voronoi diagram on the points of P, 

Z?(x) := |T/p,(x)| = |{yGr : |y - x| < |y-v,| , Vi}| 

= area of cell for x in Voronoi diagram on the points of P'. 
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By a slight abuse of notation, define 

0(x)nC, ~\Vp,{^)nVp{yr,)\. 
The notation is shown in Figure [S] The Sibson coordinates are defined to be 

^ , , D(x) n . , , T / X -Dfx) n c, 

'^n^) •= Kr~\ — equivalently, Aj(x) = —j^ . 

o(x) Ef=iO.(x)nc, 

Milbradt and Pick [5T] modify the definition of the Sibson functions for polytopes so that the coordinates of a 
point on an edge or facet of the polytope are dependent only on the Sibson functions associated to the boundary 
vertices of that edge or facet. This ensures C*^ continuity of the functions across adjacent mesh elements. 

Moreover, it has been shown that the Sibson functions are C°° on the polygon except at the vertices where 
they are C" and on circumcircles of Delaunay triangles where they are [22lll2j . Since the finite set of vertices are 
the only points at which the function is not C^, we conclude that G H^{K) where K is the domain mesh. This is 
the typical continuity required for finite element applications with nodal interpolation functions and makes them fit 
for use in the dual Whitney functions we define next. 

Definition 2 The dual Whitney function W^„3-k associated to the fc-dimensional element -ka^~^ in a 3D dual 
mesh is defined as follows. 

— Dual Vertices. The function associated to a dual vertex -ka^ := is the Sibson coordinate for the vertex, i.e. 

W,,3 := A, 

— Dual Edges. The function associated to an oriented dual edge *cr^ := [vj, Vj] is the vector-valued function 

W^^2 := A,VAj - Aj VA, 

An example is shown in Figure |4l 

— Dual Faces. Consider a dual face ★cr^ with m vertices {vq, . . . , Vm_i}. Partition the face canonically into 
triangles by adding a vertex c at the centroid of the face vertices and adding the edges [c, Vj]. Define 2-simplices 
Ti := [c, Vj, Vj^i], indices taken mod m. Define 3-simplices by connecting the to the endpoint of inside the 
polyhedron. Define 



^ — ' I ★ cr-'- 



1=0 



where is the characteristic function on r.^ (1 on x^, otherwise) and 

Wr, — 2 (AcVA; X VA,+i - A.VAe X VA,+i + A,+iVAc x VA,) . 

Note that Wt; is the Whitney 2-form associated to face of a tetrahedron (see (|26|) in Appendix A) and that 
these tetrahedra partition the entire polyhedra. An example is shown in Figure [S] 

Dual Cells. The scalar-valued function associated to a dual cell ★a" is a constant function on the cell: 



l/|*cr°l on *o-0 

otherwise 



Since the dual Whitney functions use a generalization of barycentric coordinates, it can be shown that they have 
the standard continuity across faces, e.g. tangential continuity for W^^r^ ^md normal continuity for W^^^ji . This means 
the image of Xg is in , the image of Ii is in _ff(curl), and so forth (see Figure[2]). A proof of this and other properties 
of I/j appears in [T3] . We are also developing a higher order version of these operators [IS] . 

Using dual Whitney functions, we define a novel dual discrete Hodge star by 

mk'"''r% (w,,.,w,,.). (7) 

The inner product here is the standard integration of scalar or vector valued functions over the dual domain ★/sT. For 
instance, in the case fc = 3, the definition yields 



^ • J*K 



The formulation for other k values will similarly involve integrals of the Aj functions. 
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Fig. 4 The dual Whitney function associated to the lower right edge of a pentagon is shown on the left. The magnified portion 
shows the vector field in the neighborhood of this edge. The gradients were approximated in Matlab using a simple 2-point 
difference rule on a regular grid laid over the pentagon. 




Fig. 5 Sample computation of a dual Whitney function associated to a dual face *cr^ with vertices v^. By adding the ccntroid 
c, we have a canonical decomposition of into triangles t^. A weighted sum of the primal Whitney function associated with 
each Ti is constructed to define the function for the face. As shown on the right, each t^, e.g. the shaded triangle, forms a 
tetrahedron by connecting its vertices to the vertex of cr^ interior to the polyhedron. Note that in general c need not be the 
same as cr^ D *cr^ . 



Lemma 1 (M^ is sparse. 



Proof Observe that W^^^k has localized support by construction. Entry ij of (M^"'^')""'^ will be non-zero only if ★erf 
and -ka^ are adjacent. Thus each row of the matrix will have at most as many non-zero entries as -ka^ has adjacent 
n — k cells, meaning the matrix is sparse. 



Lemma [T] does not hold if is replaced by M^'*'* as these sparse matrices typically have dense inverses. 

Note that (M^'°^)~^ is trivially sparse since it is diagonal, however, it can only be employed when the meshes are 
orthogonal. 
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3.1 Local Structure of Discrete Hodge Stars 

The continuous Hodge star * is a local operator meaning its effect on a differential form evaluated at a particular 
point on a manifold depends only on the geometry of a local neighborhood of the point. The discrete Hodge star 
is thus required to be a local operator as well meaning the evaluation of M;; on a basis cochain (1 on cr*^ and 
otherwise) should involve values on only a few simplices adjacent to erf. In the language of matrix theory, this 
requirement says M^^ should be sparse. 

We now give a more specific characterization of the sparsity structure of M^'"* and (M^"°')^^. The intuition 
for these results is demonstrated by Figure [6] 




Fig. 6 The various discrete Hodge stars depend on different aspects of mesh geometry as shown in this 2D examples. The 
diagonal Hodge star (left) computes ratios of sizes of primal-dual element pairs. The Whitney Hodge star (middle) has entries 
of Whitney functions integrated against each other. The support of a particular W^j-i function is shown in grey; the integral of 
its projection to the bold edge has value 1. The Dual Hodge star (right) that we propose has entries of dual Whitney functions 
integrated against each other. The support of a particular W^^i is shown in blue; the integral of its projection to the bold dual 
edge has value 1. 



Lemma 2 Entry ij in Mj^'"* is non-zero only if there exists cr" £ K such that a" has at least one vertex from af 
and one vertex from aj . 

Proof Computing entry ij in M^''** involves [J] Prop. 9.6] summing terms of the form 

(^l^\iX2^ det(vi^ Wj) (8) 

where Ai,A2 are barycentric functions associated to vi £ af, V2 G (Jj, respectively; 7 is a list of k vertices from af 
not including vi; J is a list of k vertices from aj not including V2; and V/, Wj are n x k matrices. The pth column 
of Vi is the vector VAp where Ap is the barycentric function associated to the pth entry in I. The qih column of Wj 
is the vector VAg where Xq is the barycentric function associated to the gth entry in J. 

Observe that the support of the barycentric function associated to vertex v is contained within the n-simplices 
touching V. Thus, if there is no cr" with at least one vertex from af and one vertex from aj, the Ai and A2 appearing 
in (jS]) will always have disjoint support, making the entry zero. 

Using the same kind of reasoning, we have a similar result for our dual discrete Hodge star. 

Lemma 3 Entry ij m (M^""')""'^ is non-zero only if there exists ★ct'' £ *A' such that ★cr" has at least one vertex 
from -ka^ and one vertex from -kaj . 

The number of fc-simplices in an n-simplex is (^^j^) which gives the following corollary. 

Corollary 1 Let A{a^) denote the number of n-simplices m K incident on at least one vertex from . Then the 
number of non-zero entries m row i o/M^''** or row i o/(M^"°')~-^ is at most (fc+i)^(o-.f)- 

The bound can be sharpened for particular choices of n and k or if additional assumptions are made about K. 
As stated, however, the corollary provides a simple means for evaluating the computational expense of a particular 
discretization scheme as we will discuss in Section 
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3.2 Numerical Stability 

To maintain the numerical stability of a DEC-based method, the discrete Hodge star matrix should have a bounded 
condition number. Put diflterently, the entries of the matrix should be roughly the same order of magnitude. This 
requirement is frequently considered from the context of numerical analysis but is often absent from the literature on 
discrete operators. 




Fig. 7 Examples illustrating how the measure of a primal simplex a'' (black) and its dual *a'' (red) need not be the same 
order of magnitude, (a) In this 2D example, the ratio | *(7^|/|(t^| can be made arbitrarily small by increasing the length of . 
(b) The ratio | *a'^\/\a''-\ can bo made arbitrarily large by decreasing the length of cr^ . (c) The ratio | *a \/\a \ can be made 
arbitrarily large by decreasing the area of . Thus, a discrete Hodge star involving terms of the form | * (t*|/|(t*| may have a 
bad condition number unless primal and dual mesh quality is controlled. 

The common thread in the geometrically-defined discrete Hodge stars such as Mf is a measurement of the 
size of dual cells i.e. | * cr'^ | . This suggests that geometric criteria on primal elements alone will not be sufficient to 
control the condition number of the discrete Hodge star matrix. In particular, since ratios of primal to dual cells are 
computed, the following criteria must be satisfied: 

Nl. Primal simplices satisfy geometric quality measures. 
N2. Dual cells satisfy geometric quality measures. 
N3. The value of | * a"'^|/|a"'^| is bounded above and below. 

N4. The primal and dual meshes do not have large gradation of elements, i.e. miuj jcr^j and maxj \a^\ are the same 
order of magnitude and min^ | *cr,f | and max^ |cr,f | are the same order of magnitude. 

Conditions Nl and N2 are required for discretization stability. Aspect ratio is often used as a geometric quality 
measure for tetrahedra. Conditions N3 and N4 are based on our analysis above. Condition N4 in particular shows 
that these discrete Hodge stars are not fit for use on meshes tailored to multi-resolution situations where gradation 
is necessary to achieve reasonable computation times. Examples are shown in Figures [7] and [H] 




Fig. 8 Graded meshes also present a problem for discrete Hodge stars involving primal-dual size ratios. The primal mesh shown 
here induces a wide variation in values of 

I *o-*|/|o-*| for k = 0,1,2. This can cause ill-conditioned Mj. matrices, resulting in numerical instability. 

For Mf^'\the size of the matrix entries are controlled by the size of the inner products of Whitney basis forms. 
The integrals in ^ are on the order of the size of |(Tj,|, meaning again that a large gradation in primal mesh element 
size could produce large condition numbers. Since M^'"* does not depend on the size of dual mesh elements, however, 
its condition number is more stable against violations of conditions N2 and N3. Analogously, the condition number 
of (M^"°')~^ is more stable against violations of conditions Nl and N3. Our conclusions are summarized below. 
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— Conditions N1-N4 are necessary to ensure '"^ has a good condition number. 

— Conditions Nl and N4 are necessary to ensure M^^^^ has a good condition number. 

— Conditions N2 and N4 are necessary to ensure (M^""')^^ has a good condition number. 



3.3 Improved Condition Numbers with (Mf ""')"^ 



To provide concrete evidence for our numerical stability claims, we present a simple example in 2D showing how 
M, " and Ml can have condition numbers an order of magnitude worse than (Mf"°')~l on the same mesh. 
This serves as a proof of concept that the DEC-based dual formulation of a problem can provide practical advantages 
in cases of difficult mesh geometry. 




In the 2D mesh shown in Figure[51 the labeled vertices of the primal mesh have coordinates vi = (0, 0), V2 — (0, 1), 
V3 — (P, and V4 — (— P, where P is a free parameter we can adjust to modify the geometry. The remaining 
vertices are chosen so that they form equilateral triangles with edges (713, 0-23, cri4, and (724, as shown. The orthogonal, 
circumcenter-based dual mesh is shown in red. 

Without loss of generality, fix any ordering on the mesh edges, beginning with 



{0-12, 0-13, 0-14,0-23, 0-24, •■ •}• 



(9) 



We first calculate the upper left 5x5 block of each matrix, yielding the matrix values assigned to all possible 
interactions between pairs of these first five edges. Using the circumcentric dual mesh and definition Q , we compute 



/ 4P' 



^Diag 



4P 


- 







p 





p 





£1 





£1 



\ 



(10) 
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Since 



fDiag 



is diagonal, its condition number is the ratio of its largest diagonal entry 



where q^^ + -1 
to its smallest. The uncomputed diagonal entries will be very close to q meaning the condition number can be 
approximated as 



cond 



10 



Using the Whitney interpolant for edges (see ()25|) in Appendix A) and the definition of M]^'"* given in ([S]), we can 
also compute 



^7050 
/3 7 5 
/? 5 7 
^ (5 7 



V 



(11) 



7 



where a 



P = 



suggested by is an artifice of our ordering of the edges as stated in ([9]). However, the remaining diagonal entries 
of M]^'"* are all close to 7, the entire matrix is symmetric, and the remaining non-zero off-diagonal terms are roughly 
the same size. Thus, the eigenvalues of the 5x5 matrix shown in (jlip allow us to approximate the condition number 
of M];^'''*. Using Mathematica, we find analytical expressions for the max and min eigenvalues of the 5x5 matrix 
and take their ratio to approximate 



12P^+20\/3P+21 



and & — 



4P^-5 



Note that some of the structure of . 



vWhit 



cond 



24P2 + 5\/3P -f a/288P4 - 120\/3P3 + 3p2 _^ 9 ^_ 3 



eO(P) 



10\/3P + 18 

Finally, we compute (Mf using the barycentric dual mesh and definition (O, yielding 

C 61 fi: C •■ 



C K e C 

C ? e fi: 



C c ^ 



V 



(12) 
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where i9 = (^ry^^i^ , 7?^^i J , C = (^t?^,,!^ , J7^<,i J , e = (^r?^^i^ , ry^^i J , k = (^77^^!^ , ry^^i J and ^ = (^77^^!^ , ry^^i J . Note 
that analytical expressions of these inner products are not feasible due to the need to compute areas of intersection of 
irregular polygons in the definition of the A functions. Instead, using Matlab, we create a simple grid-based quadrature 

method to estimate the entries of ^Mf ""'^ for various values of P. As with M]^'"* , we then estimate the condition 
number of the entire matrix by the ratio of the max and min eigenvalues of the 5x5 matrix given in (|12p . 

The cases P = 2, 5, and 10 were tested. The integral required to compute ^ has support outside of the portion 
of the dual mesh shown in Figure |9] We thus set ^ to be the same as C,, since both are inner products associated to 
adjacent edges in the dual mesh. The computed values of k were very small, as expected; we found that setting k to 
zero did not affect the condition number estimate. Our results are summarized in Table [T] 



p 


cond (Mf 


cond {mY^^'^) 


cond((Mf"°') ^) 


2 


6.3 


3.2 


1.5 


5 


17.2 


9.9 


1.3 


10 


34.6 


21.6 


1.4 



Table 1 Comparison of condition numbers of different discrete Hodge stars for various values of P. 
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primal: 



dual: 



Fig. 10 Portion of a generic DEC-deRham diagram (cf. Figure [2]l showing the natural duality between the variables and 
operators of systems II13I I and II14II . Discretizations of the variables are written in place of the primal or dual cochain spaces to 
which they belong. 



Our numerical experiments thus provide evidence for the claim 

cond(^(Mf""')-^^ G 0(1). 

The above example confirms that while our dual discrete Hodge star has an analogous definition to the primal 
discrete Hodge star, its condition number is indeed controlled by the geometric properties of the dual mesh elements, 
not those of the primal mesh elements. This fact is especially useful for problems on tetrahedral meshes where slivers 
(narrow, nearly planar tetrahedra) frequently occur and are difficult to remove. 



4 Applications 

The dual interpolation functions In-k we defined in (O and the dual discrete Hodge star we defined in (O are new 
tools for designing stable finite element methods. We start by explaining the generic methodology of our approach 
and then apply it to two sample finite element problems from the literature: magnetostatics and Darcy flow. 



4.1 Generic methodology 

The Discrete Exterior Calculus approach to discretizing a PDE is as follows: 

I. Translate the continuous PDE problem into the language of exterior calculus. 
II. Linearize the problem, possibly by introducing an intermediary variable (i.e. a mixed method). 
HI. Discretize the fc-forms into fc-cochains and the operators d and * into D and M matrices. 
IV. Solve a linear system constructed from the discrete equations. 

Our methodology focuses on step III and exposes how there are often many natural choices for discretization in line 
with DEC theory. Consider the case where we are given a PDE in terms of a variable u that is treated as a fc-form 
in the continuous setting. Suppose that a mixed method is possible in which the intermediary variable v should be 
interpreted as an n — A; — 1 form. In this case, the typical mixed linear system is 



V \ / F 



(13) 



J \v J \G 

where u G C*^, v G C" ^ ^ axe the discretized variables and F €C"^ ^ , G € C^^^ are the discretized load data. 

The simple idea at the heart of our technique is to swap the type of dicretization (primal or dual) of each variable 
and then infer the rest of the system from DEC theory. Note that the cochain order of each variable should not 
change, only the mesh on which it is discretized. Hence, the dual formulation of system (|13p is 



(14) 



where now u G C^, v G ^"^'^'^1 are the discretized variables and F G C"^'^, G G c''^^ are the discretized load 
data. We show in Figure [10] how these two discretizaions flt into a generic DEC-deRham diagram in a natural and 
complementary fashion. 

Additional equivalent systems can be derived by using proxy variables in clever ways, e.g. solving for some 
z G C'^"^ such that x is deflned uniquely by x = Dfe_iZ. These systems are easiest to understand via the speciflc 
examples we now examine. 
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4.2 Magnetostatics 

The magnetostatics problem is characterized by Gauss's law for magnetism, Ampere's law, and a constitutive rela- 
tionship, respectively, 

div 6 = 0, *b — h, curl h — j. (15) 

Here, j is a given current density and b and h both represent the magnetic field. It is assumed that the domain O 
is contractible with boundary F written as a disjoint union U F^ such that n ■ b — on F'^ and n x ft = on f'^. 

A DEC-based treatment of the problem reveals canonical and symmetrical ways to put this into a mixed formu- 
lation linear system, depending on whether b is discretized as a primal or dual cochain. If we discretize 6 as a primal 
2-cochain B G and /i as a dual 1-cochain H £ C , equations p5|l become 

D2B = 0, M2B = H, bJb = J. 
This allows for two possible mixed systems. The first is 



I2 W B \ _ / -Ho 



p \ 



(16) 



In this system, Hq £ is any dual 1-cochain satisfying D^Hq = J and H is defined by H := Hq + D^^P. Thus 



l)f H — I$f{liQ + D^^p) = J is assured. 
The second mixed system is 



h\ / 



(17) 



In this system, B is defined by B ;= DiA, so that D2B = D2D1A — 0. For a fixed j, systems (|16p and p7|l result in 
the same solution pair (b, ff) and were shown by Bossavit [7] to converge to the solution pair (b, h) to (|15p as the size 
of mesh elements goes to zero. 

2 

We now consider a novel dual discretization approach by treating b as a dual 2-cochain BSC and /i as a primal 
1-cochain l-l £ C^. The continuous problem H15|) is now discretized by 

DoB = 0, B=MiH, ]D)iH = J. 

The first mixed system of this dual formulation is 

Do 



(18) 



In this system, Hq £ is any primal 1-cochain satisfying DiHq = J and H is defined by H := Mj^ ^B. Thus DiH 
Di(ho -I- DgP) = J is assured. The last system is 



(19) 



where B is defined by B := D^A so that D^B = DqBJa = 0. For a fixed j, systems ((TSJ and (fl^ will result in the 
same solution pair (b, h). In a future work, we will show that these systems also converge to the solution pair {b,h) 
to (|15|l as the size of mesh elements goes to zero. Taking that for granted, we state the advantages of having all four 
systems (|16|l . (|17|l . (|18|l . and (|19|l available for implementation. 

First, observe that systems (|16p and (I17|l make use of the M2 matrix and its inverse while psp and (|19p use 
the Ml matrix. If the diagonal Hodge star is used, then M2 requires good ratios between the size of primal faces 
and their dual edges while Mi requires good ratios between the size of primal edges and their dual faces. Thus, on 
unstructured meshes, one system may break numerically on a mesh that is acceptable for another system. 

Second, if the Whitney Hodge star is used, M^"'^ may be a full rank matrix, making systems (|17p and (|18p less 
attractive numerically. By constructing the dual discrete Hodge stars as proposed in this paper, these systems become 
sparse again by Lemma [T] and thus are available as a practical alternative. 

Third, having four systems available for the same problem allows for rigorous error- checking and cross-confirmation 
of results. This is particularly valuable when physical experimental confirmation of the results is impossible or 
expensive. 
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4.3 Darcy Flow 

The Darcy flow problem in is 

f / + |vp = in 

< div / = in (20) 
l_ / ■ n = ?/) on dfl, 

where fc and are physical constants, / is volumetric flux and p is pressure. It is assumed that there is no external 
body force, the boundary F := dfl is piecewise smooth, and the compatibility condition J^^ = J^^^ -i/jdF is 
satisfled. Without loss of generality, take ^ = k. 

First consider discretizing / as a a primal 2-cochain F G C and p as a dual 0-cochain P G C , yielding the 
discretized equations 



Hirani el al. [20] used this approach to derive the linear system 

/ n \ 

(21) 



/ V p / V* 



We present an alternative formulation using the same discretization, inspired by the magnetostatics systems (|17p and 
(|19|l . Let Fq G be a primal 2-cochain satisfying D2F0 = ^. The system is 



vf-i Di\ /q\ ^ /-Fo 



(22) 



Here, P is a solution to DI'p = — Q. The existence of P is guaranteed by the exactness of the dual cochain sequence 
at C and uniqueness of P is determined by initial conditions or boundary data. The flux cochain F is defined to be 



n^^Q so that D2F = D2(M^^Q) = D2(fo + DiG) 



_ 2 

We now present the dual formulations derived by treating / as a dual 2-cochain F G C and p as a primal 0-cochain 



p G C*'. The discretized equations are now 



The first system of this formulation is 



lIj"^F-HDoP = 0, ^}l¥ = $. 



The second system is 



(23) 



(24) 



where Fq is a solution to Dq'^o = ^ a^nd F is defined to be MiQ, analogous to system H22[) . Thus, taking Dj^ of both 
sides of the top equation of (|24p yields DJ^f = Further, the bottom equation of (|24p yields DiQ = which, by the 
exactness property of the primal cochain sequence implies that there exists a solution p to DqP = ~Q- 

We now have four mixed systems, (|21l) - H24p . discretizing the Darcy flow equations (|20p . three of which had not 
be considered by Hirani et al. [20] . This plethora of equivalent systems offers the same advantages as those discussed 
at the end of the magnetostatics example from Section [4.21 



5 Conclusion 

In this work we have augmented the theories of Discrete Exterior Calculus and mixed methods by introducing two 
novel tools: Whitney-like interpolation functions defined on dual domain meshes and a sparse inverse discrete Hodge 
star. We have shown the tools to have natural, straightforward definitions and clear geometric interpretations. We 
have used them to derive previously unexamined numerical stability criteria relating to the condition number of the 
discrete Hodge star used in the method, based on the geometry of the dual mesh cells. Further, we have demonstrated 
in both general and specific contexts how these tools can be used to develop alternative discretizations of PDEs with 
sparse, well-conditioned matrices. The techniques we have described provide a valuable methodology for researchers to 
revisit their current finite element formulations and confirm or improve their results with new discretization methods. 
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A Whitney Functions for Primal Meshes 

Whitney fc-forms arc picccwisc linear functions on a primal mesh, one for each fc-simplex in the mesh. 

— Primal Vertices. The Whitney O-form associated to a vertex cr" := is denoted 

where is the barycentric function for the vertex. More precisely, is defined by the condition of being linear on every 
simplex of the mesh, subject to the constraints Ai(vj) = Sij. 

— Primal Edges. The Whitney 1-form associated to an oriented edge := [vi,Vj] is the vector-valued function 

W^i := A,VAj -A. VAi. (25) 

— Primal Faces. The Whitney 2-form associated to an oriented face := [vi,Vj, v^] is the vector-valued function 

W^2 := 2 (AiVAj X VAfc -I- Aj VAfc X VA^ + AfeVA, x VA^) (26) 

— Primal Tetrahedra|3 The Whitney 3-form associated to an oriented tetrahedron is its characteristic function, scaled 
by the reciprocal of the volume a"^ . 

l/|o-3| on 0-3 

otherwise 



B GeneraHzed Barycentric Functions 

Let 7" be a top-dimensional cell of the dual mesh (i.e. a polygon in 2D or a polyhedron in 3D) with vertices vi, . . . , vjy- A set 
of functions A; : T — > K, j = 1, . . . , Af are called barycentric coordinates on T if they satisfy two properties. 

Bl. Non-negative: A^ > 0. 

B2. Linear Completeness: For any linear function L : T — > M, 

AT 

L = ^L(v,)Ai. 

i = l 

A set of barycentric coordinates {A^} also satisfies these additional familiar properties: 

AT 

B3. Partition of unity: A^ = 1. 

1=1 
JV 

B4. Linear precision: Vj Xj (x) = x. 

i = l 

B5. Interpolation: Ai(vj) = Sij. 

A proof that properties B3-B5 arc implied by B1-B2 in the 2D case can be found in our paper 1161 . The 3D case is similar. 

Three major approaches to defining generalized barycentric functions on 2D polygons have emerged in the literature. The 
Wachspress functions I24II13I are rational functions constructed explicitly based on the areas of certain triangles within T. The 
Sibson functions 1221 . also called the natural neighbor or natural element coordinates 1231 . are also constructed explicitly, but 
instead use the areas of Voronoi regions associated with the vertices of T. The Harmonic functions I26l l8l are defined as the 
solution to Laplace's equation over T with certain piecewise linear boundary data. 

We have shown in 1161 that any of these functions suffice to give the optimal interpolation estimate for the lowest order 
case in 2D, assuming some basic geometric quality criteria on the dual mesh elements. For this paper, we have employed only 
the Sibson coordinates as they generalize easily to 3D, are reasonable to implement, and are more stable against bad geometry 
than the Wachspress functions. A formal proof of their convergence properties in 3D will be the focus of a future work. 

Acknowledgments We are grateful to Alexander Rand for his help in implementing the Sibson coordinates. This research was 
supported in part by NIH contracts R01-EB00487, R01-GM074258, and a grant from the UT-Portugal CoLab project. 



^ Note that the W^3 definition has been simplified from a more general definition of Whitney forms 1271 using the geometric 
identity 

VA, ■ (VAj X VAfc) = ±-^37 

where the right side has sign —1 if an odd index was omitted from the scalar triple product and -f 1 otherwise. This reduces the 
sum in the general formula to (l/lc'* |) A;, which is simply l/lc^l due to the partition of unity formed by the barycentric 
functions. 
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